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Abstract 

In this paper we report on 2D numerical simulations concerning linear and nonlin- 
ear evolution of surface-tension-driven instability in two-fluid systems heated from 
below using classical and phase-field models. In the phase-field formalism, one intro- 
duces an order parameter called phase-field function to characterize thermodynam- 
ically the phases. All the system parameters are assumed to vary continuously from 
one fluid bulk to another (as linear functions of the phase- field) . The Navier-Stokes 
equation (with some extra terms) and the heat equation are written only once for 
the whole system. The evolution of the phase-field is described by the Cahn-Hilliard 
equation. In the sharp-interface limit the results found by the phase-field formalism 
recover the results given by the classical formulation. 
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1 Introduction 



In systems with multi-material interfaces the classical methods may admit 
some non-physical oscillations near interfaces or meet with difficulties in the 
problems with large deformations. Therefore recently R. Fedkiw et ai devel- 
oped a new numerical method— the so-called "ghost-fluid" method— for de- 
scribing interfaces in multi-material flows [Caiden et a/., 2001], [Fedkiw et 
al.^ 1999], [Fedkiw & Liu, 2002]. In this method a so-called level set func- 
tion (which satisfles the level set equation) is implemented to keep track of 
the interface. The zero level marks the location of the interface, the positive 
values correspond to one fluid and the negative values to the other. In addi- 
tion they captured the appropriate interface conditions by deflning a ghost 
fluid (for each fluid) which has the same pressure and velocity as the real 
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Fig. 1. Sketch of the system: two superposed immiscible fluid layers. The system 
is heated from below and the temperatures on top and bottom are maintained 
constant. 

fluid, but the entropy of the other real fluid. Since the ghost fluids have the 
same entropy as the real fluid which is not replaced, a one phase problem is 
solved for each fluid. Another continuum model adequate to multi-systems 
in which the interface location cannot be explicitly tracked is based on the 
phase-field method. Proposed for the first time by Langer [1986], this model 
was subsequently developed for studying solidification phenomena [Anderson 
et al, 2000], [Langer, 1986], [Nestler et al, 2000], [Wang et al, 1993], den- 
dritic crystal growth [Braun & Murray, 1997], [Tong et a/., 2001] or dynamic 
fractures [Karma et a/., 2001]. In this type of model, in place of the set level 
function an auxiliary variable— the phase-field (p —is appended to the usual 
set of thermodynamic variables in order to provide an explicit indication of 
the thermodynamic phase in each point of the system. This order parameter 
(p can be also regarded as a fluid concentration in a binary mixture of two 
fluids. The phase-field takes distinct values in each bulk phase and undergoes 
a rapid but smooth variation in the interfacial region, (p is governed by a par- 
tial differential equation that guarantees the realistic interface conditions in 
the limit of a suitable thin interface. Therefore, the problem is treated like 
an entire (one phase problem), exactly as in the ghost-fluid model. Unlike the 
ghost fluid model, in the phase-field formalism the interface is diffuse with a 
finite thickness over which the fluid properties vary smoothly from one region 
to the other. Consequently a new phenomenon is taken into discussion: the 
diffusive transport in the interfacial region between the two phases. 
Hence, the continuum models are suitable for describing multi-layer systems, 
because they reduce the system of equations (they don't need different equa- 
tions for each medium) and avoid the interface conditions. 

Due to these two essential advantages, we will adjust in this paper the phase- 
field model to another problem with material interface: the two-dimensional 
Marangoni convection (MC) induced by the temperature gradient in a two- 
layer system of immiscible and incompressible fluids (see Fig. 1). As it is 
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Fig. 2. Phase-field distribution versus z for the stationary state. At = (fluid 

1 boundary) the phase- field function takes the value (p = —1 and at z = 2 (fluid 

2 boundary) the value (f = +1. The diffuse interface between the two different 
immiscible fluids is situated around z = 1. 



well known, the equation which describes the bifurcation from the conductive 
regime to the convective regime for the system illustrated in Figure 1 can be 
casted into the following normal form [Colinet et al, 2001], [Manneville, 1990], 
[Newell k Whitehead, 1969]: 



— ^{(5 + i^)A-a\A\'A 



where A denotes the amplitude measuring the " intensity" of convection and [3 
is the control parameter, connected with the temperature gradient. (If Re{a) 
is negative, higher order terms have to be included in the above equation in 
order to ensure the global stability.) When 7 = the bifurcation from the con- 
ductive regime (/? < 0) to the convective regime (/? > 0) is realised through a 
monotonic instability and when 7 7^ through an oscillatory instability. For 
the last situation, the Poincare map [dA/ dt^ A) consists in a limit cycle, and 
this kind of bifurcations is known in literature as Hopf bifurcation [Manneville, 
1990]. 

In our paper we will report on numerical simulations concerning linear and 
fully nonlinear evolutions of surface-tension driven instability (or MC with 
short wavelength) , corresponding to the two different bifurcations types men- 
tioned above. After the introducing of the method, we shall discuss a silicon- 
oil— air system heated from below, a situation for which the short-wave in- 
stability appears as monotonic mode. Then we particularize the results for 
a water— n-octane system, heated also from below. Here the surface-tension- 
driven instability appears as oscillatory mode via a Hopf bifurcation. The 
results given by the phase-field formalism will be compared with those given 
by the classical model, assuming flat interface and two sets of equations for 
both fluids coupled by interface conditions. 



3 



1 r 



-1 



-2 



-1 



2 



Fig. 3. Free-energy density representation versus phase-field function ip. We have 
chosen f{(p) as a symmetric " double- well" potential with two local minima: one at 
= corresponding to the fluid 1, and the second at = +1, corresponding to 
the fluid 2. 

2 The model and basic equations 

The phase- field is assumed to be = — 1 at the lower boundary of fluid 1 
{z = 0) and (/9 = +1 at the upper boundary of fluid 2 (z = 2) while the 
diffuse interface between the two different immiscible fluids is situated around 
z = 1 for the simulations examined in the present paper (Fig. 2). The system is 
bounded in the vertical direction by two solid, perfectly heat conducting walls 
with fixed temperatures T{z = 0) = and T{z = 2) = Tf. In the phase-field 
formalism the Navier-Stokes (NS) equation must contain some non-classical 
phase-field terms. In the following we shall derive these supplementary terms 
by minimizing the free-energy functional for our system [Anderson et al.^ 1998], 
[Jasnow & Vinals, 1996], [Langer, 1986]: 



The first term from the above functional represents the free-energy density for 
an homogeneous system corresponding to regions far from interphase bound- 
aries. The latter one characterizes the energy concentrated in the interfacial 
region. f{Lp) is considered to be given by the expression f{Lp) = — ^) 
(C— positive constant) [Borcia & Bestehorn, 2003] which is a symmetrical 
curve with two local minima: one corresponding to (/^ = — 1, namely for the 
fluid 1 bulk, and another one to 99 = +1, for fluid 2 (see Fig. 3). The coefficient 
/C in (1) is connected to the surface tension coefficient, because the second 
term from the functional (1) is associated with variations of the density (and 
consequently of the phase-field) and contributes to the free-energy excess of 
the interface [Jasnow & Vinals, 1996], [Rowfinson & Widom, 1982]: 
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In many previous works /C is assumed to be constant [Anderson et a/., 1998], 
[Jasnow & Vinals, 1996], [Langer, 1986]. But the short-wave instabihty is 
driven by the gradients of the surface-tension. To include the thermo-capiUary 
effects it is necessary to consider /C weakly dependent on temperature: 

K^K^-KtT {JCt>0). (3) 



The Lagrangian energy density jC of our system which is assumed to have 
an explicit dependence on the scalar its gradient V(p and on the gradient 
energy coefficient /C(T) 

c {if, Vvp, /c(r)) = /(vp) + ^{Viff (4) 



must satisfy the Euler-Lagrange equation: 



dC ' d ( dC \ 

dip Udx,\d{d,^)) • 

Substitution of relation (4) in (5) gives us the following relation (Noether's 
theorem) : 
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d f dC dip dCdIC, 



dxi \d{diip)dxj J dICdx, 



= 0. (6) 



Assuming K. = const., the stress tensor is defined as: 
S = ICVif ®V(p-C I, 

or, in components, 

- - - rs 

dxdx- 

and (6) can be considered as a conservation law: 

V • H = 0. (7) 



Under these conditions, a conservative term — V • 5 appears in the Navier- 
Stokes equation. This term describes the contribution of capillary forces [An- 
derson et al, 2000]: 
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dv 

p— = V{-p + C)-V- {JCVip ® V(p) 

(JjL 

+V • {r]\/v) + pgz. (8) 

(Both fluids are assumed incompressible.) But if /C 7^ const. ^ the conservative 
law (7) is not satisfled. Hence we must consider also the last term of relation 
(6). Therefore a new force-component has to be included: 



p_ = V (-J9 + £) - V • (/CV^ ® V<^) 

+V • {rjVv) + pgz + \lCTVT{V^f. (9) 

Li 

This last term from (9) represents the "Marangoni force" which is the trigger 
for the short-wavelength instability, because integrating equation (9) through 
the interface, this new-force component assures the fulfllment of the inter- 
face condition for tangential stresses in the limit of sharp and rigid interfaces 
[Borcia & Bestehorn, 2003]: 

^L(2) -^L(i) = 



Wxz = ^ + ^) the viscous stress tensor for incompressible fluids). 
Equation (9) together with the energy equation 

dT 

pc— = v-{Kvn (10) 



(p is the fluid density, c— the heat capacity and tv—the thermal conductivity) 
and the Cahn-Hilliard (CH) equation, which serves for mass conservation in 
a system with diffuse interface 

^^-MoAllCAip-^] (11) 



(Mo— the phenomenological mobility), describe both Marangoni instabilities 
without to impose supplementary interface conditions. In the classical models 
the basic equations are written two times: once for the fluid 1 and the second 
time for the fluid 2 and matched by interface conditions. In our phase-field 
formalism the NS equation (9), the heat equation and the CH equation are 
valid for both fluids while all the system parameters are assumed to be linear 
functions of the phase-field: 

Pi + P2 Pl- P2 Vl +V2 Vl- V2 

— :^ :^ — V = — ^ ^ — ^ 
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Ci + C2 Ci - C2 /^i + 1^1 - 1^2 
C — (Z?, Hi = (ZP. 

2 2 ^' 2 2 ^ 



In the above relations the index "1" describes the fluid parameters at z = 0, 
while the index "2" describes the fluid parameters at 2; = 2. After adimen- 
sionalization of (9)-(ll) 

T , dl , _ xi - , T — Tt 

r = dir. t = — t , = — , T = — — -, 

Xi ^1 76 -T,' 



(Xi— thermal diffusivity at the lower boundary of fluid 1, di— the fluid 1 depth) 

the following non-dimensional parameters appear: 

Pr — — — Prandtl number of the fluid 1, 

Ca = ^y^di _eapiUary number, 

M = ^TvTO m-Tt)di _Marangoni number, 

G = Galileo number, 

L — ^y^— the width of interface 

M'^ = the phenomenological mobility. 

Assuming both fluids being incompressible, the velocity fleld, v{x^ z, t) can be 
replaced by the stream- function, ijj{x^z^t) (the primes are dropped): 



with i and /c the unit vectors in x and z direction, respectively. 



3 Purpose 



For small capillary numbers Ca and Galileo numbers G, that is, for very thin 
layers, the interface between the two immiscible fluids is deformable. Thereby 
the long- wave instability induced by surface deflections appears even for small 
Marangoni numbers. For large capillary and Galileo numbers the interface 
becomes quasi-non-deformable and in the system develops only one instability: 
the short-wave instability. The goal of our paper is to analyze this instability 
in the frame of the phase-field model for the following situations: 

1. the interface is perfectly rigid but diffuse (model 1); in this case we drop CH 
equation (11) and we assume for the phase-field in (9) a particular variation 
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Fig. 4. Dependencies of critical Marangoni number Mcr and critical wavenumber 
kcr on the interface thickness, for MC driven by surface tension gradient. The plots 
correspond to two different models: model 1 assumes a perfect rigid liquid-gas in- 
terface, while in model 2 the interface is quasi non-deformable. For model 2 one 
considers: = 2 • 10^ G = 3 • 10^. Both models converge in the limit of sharp 
interfaces to classical results obtained for systems with fiat interface. 

of the form: 

ip^^\z) =tanh 



2. the interface is quasi-non-deformable; for this situation the full system of 
equations (9)- (11) is considered, but in the hmit of large values for capillary 
and Gahleo numbers (model 2). 

The results described by these two situations will be compared with the results 
given by the classical model, when the NS equation and the heat equation are 
written for both media [Engel & Swift, 2000]. At z = 1 the interface conditions 
are given by [Colinet et al, 2001], [Engel & Swift, 2000]: 

oz oz oz oz 

dz'^ PiVi dz^ dx^ ' 



where the subscripts 1(2) denotes the lower (upper) fluid, 6 represents the 
temperature perturbation, z/— the kinematic viscosity and M' is the usual def- 
inition of the Marangoni number [Engel & Swift, 2000]: 
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Fig. 5. Growth rate of surface-tension-driven-instability, Re A versus k near thresh- 
old as it results from classical and phase- field models. The three curves agree well 
near k = 2. 



M' 



(12) 



(with Bi = Kl/a—ihe Biot number, k.' = 1^2/ i^i and a = 6/2/^1)- 



4 Results 



Results concerning the linear stability analysis for the three models in bi- 
dimensional (x, z) representation are emphasized in Figs. 4 and 5 for a IQcS 
silicon-oil— air system [pi — 9A0 kg/m^^ p2 — 1.2 kg/m^^ ui — 0.1 cm^/s^ 
U2 = 0.015 cmVs, /^i = 0.142 J/msK, K2 = 0.026 J/msK, c^i = 1450 J/kgK, 
Cp2 = 1042 J/kgK [Golovin et a/., 1997]) heated from below with a — 1. For 
this stability analysis the small perturbations are assumed as 



T{x,z,t) 



( ^(o)(^)\ 

(z) 
rW(z) 



+ 



rW(z) 



exp{ikx)exp{Xt) 



with fc— wavenumber and A— the (complex) growth-rate. The derivatives with 
respect to z are expressed using a finite-difference method. Figure 4 represents 
the critical Marangoni number and the critical wavenumber versus 1/L for the 
two models 1 and 2. One can see how in the limit of sharp interfaces the results 
given by the phase- field models saturate for large values of 1/L. The critical 
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Fig. 6. Time series for the perturbation profiles in the surface-tension-driven case as 
it results from the classical formalism. The panels (a) correspond to the temperature 
perturbations and the panels (b) to the stream-function perturbations. Starting 
from an initial random pattern the perturbation evolution is represented till the 
saturation state. 

Marangoni number saturates around Mcr = 750, which means in the terms of 
the usual definition (12) a value around M^^ = 110. 



This value is in good agreement with the critical Marangoni number obtained 
by the classical model for the same lOcS^ silicon-oil— air combination with 
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Fig. 7. Same as Fig. 6, but in the phase-field description with rigid interface (model 
!)• 

a = 1: M'^^ = 93. If one compares the critical wavenumber given by the mod- 
els 1 and 2 in the sharp-interface limit with the critical wavenumber given 
by the classical model, the agreement is much better. For model 1 one ob- 
tains: kcr = 2.06, for model 2: kcr = 2.04 and from the classical model one 
obtains for the critical wavenumeber kcr = 2.003. Figure 5 shows the growth- 
rate of Marangoni instability with short-wavelength at threshold versus the 
wavenumber k as it results from classical and non-classical models. One can 
see a good concordance between the curves obtained with the classical model 
and those obtained with the phase-field model for sharp interfaces, especially 
around the threshold. For model 2 the interface between the liquid and the 



11 



C\2 I 1 

2 4 6 

X 

C\2 I 

10 



t = 2.5 



2 4 
X 

(a) 



X 

t > 45 



2 4 
X 

(b) 





2 4 6 




X 






o 


O 




2 4 6 




X 




(c) 



Fig. 8. Same as Fig. 6, but for model 2, when the hquid-gas interface is 
quasi-non-deformable. The panels (a) display the perturbation profile of the tem- 
perature, panels (b) of the phase-field and panels (c) of the stream-function. 

gas is quasi-non-deformable, but not completely rigid. For this situation, the 
curve ReX = f{k) indicates the second Marangoni instability, the long- wave 
instability, which develops around = in systems with a deformable inter- 
face. 



In order to describe the nonlinear evolution of short-wave instability for both 
diffusive models and for the classical model, the corresponding systems of 
equations were numerically solved. We used a semi-implicit scheme: implic- 
itly for the linear part and explicitly for the nonlinear one. A fast Fourier 
transform was applied in x direction and the derivatives with respect to z 
were expressed in finite-differences. To avoid the convolution sums in Fourier 
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space, the nonlinearities are calculated in real space using again a finite differ- 
ence method with grid spacing Ax and Az for their derivatives. Figures 6, 7 
and 8 display bi-dimensional (x, z) representations of the system perturbations 
for each model for the same oil— air combination and for the same s 

_ M - Mcr 

8=0.2. All the plots in the following concerning phase-field model correspond 
to large values for 1/L. Starting from an initial random distribution of very 
small amplitude, the pattern evolution is shown till saturation (moment em- 
phasized by the contour lines) . The time series were calculated on a mesh with 
256 X 40 points for the classical model, 64 x 37 points for model 1 and 64 x 51 
points for model 2. For the frames depicted in Figs. 7 and 8 (corresponding to 
the models 1 and 2) the integration time step is At = 0.1, and for those from 
Fig. 6 (classical model) At = 10~^. Choosing an iteration time step At = 0.1 
for the classical model, one obtains very quickly (after a few iteration) the 
specific pattern for short-wave instability (two-convective motions developed 
in the liquid and in the gas). The presence of a diffusive interface decelerates 
the process of pattern formation for the short-wavelength instability and also 
the process of saturation compared to the classical case where the interface is 
well defined. If additionally the diffuse interface becomes weakly deformable, 
the processes of pattern formation and saturation become longer. For the lat- 
ter situation one can also follow the time evolution for the phase-field which in 
fact describes the time evolution for the small deflections of the interface (see 
Fig. 8-b). These small deformations of the interface affect the aspects of the 
convective motions. In this case the deflections allows a momentum transfer 
from the liquid to the gas and vice-versa. Because the viscosity in the liquid 
is one order of magnitude larger than the viscosity in the gas, it follows that 
the convection in the gas becomes stronger in comparison to the convection in 
the hquid as indicated in Fig. 8-c. For all the situations represented in Figs. 
6, 7 and 8 the short-wave instability appears as monotonic instability. 



An oscillatory instability (corresponding to a Hopf bifurcation) appears in a 
water— n-octane system (pi = 998 kg/m^^ p2 = 704 kg/m^^ Ui = 0.01 cm^/s^ 
= 0.00813 cm^/s, ki = 0.6 J/msK, K2 = 0.15 J/msK, Cpi = 4182 J/kgK, 
Cp2 = 2091 J/kgK [Colinet et al^ 2001]) heated from below with a = 1. In 
this case the critical Marangoni number is: Mcr = 1.5 • 10^ and the critical 
wavenumber: kcr — 1.9. With these parameters and s = 0.2 we have illus- 
trated in Figs. 9 and 10 some snapshots till saturation state is achieved. The 
frames depicted in Fig. 9 correspond to the classical model (with a mesh of 
256 X 40 points) and those from Fig. 10 to the model 2 (with a mesh of 
32 X 51 points). In the saturated state the patterns for both models are stand- 
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Fig. 9. Time series for perturbations corresponding to short-wave oscillatory insta- 
bility as it results from the classical formulation. The left panels (a) describe the 
temperature perturbation and the right ones (b) the stream-function perturbation. 
Starting again from a random initial distribution the first two sequences depict 
the formation of short-wave instability. The last two sequences show the evolution 
during one complete period when saturation is already reached. 



ing waves. More generally, model 2 is able to reproduce the time evolution 
of the stream-function and temperature perturbations described by the clas- 
sical formalism; moreover one can follow the interface dynamics during the 
oscillation processes, a big advantage of the phase-field model. 
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Fig. 10. Same as Fig. 9, but in phase-field description with quasi non-deformable 
interface (model 2): (a) temperature, (b) phase- field and (c) stream- function. 

5 Conclusions 



We developed a phase-field model for Marangoni convection in a two-layer- 
system with diffuse interface. We reported on 2D simulations concerning surface- 
tension-driven instability via phase-field and via classical models. Linear and 
nonlinear behavior for this kind of instability is examined. In contrast to the 
classical formulation, the phase-field model treats the problem continuously, 
no interface needs to be tracked. In the sharp-interface limit our results em- 
phasize how the phase-field model applied on MC recover the results given by 
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two-layer-systems with interface conditions. Using the same phase-field for- 
malism, we plan to extend our research to the long-wave instability induced 
by surface defiections and to the infiuence of evaporation on Marangoni con- 
vection, problems which can be discussed in this approach in a more natural 
way. 
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